Imports
library(evd);
library(evdbayes);
library(coda);
library(ismev);
Lade nötiges Paket: mgcv
Lade nötiges Paket: nlme
This is mgcv 1.8-24. For overview type 'help("mgcv-package")'.
library(xts);
Lade nötiges Paket: zoo
Attache Paket: ‘zoo’
The following objects are masked from ‘package:base’:
as.Date, as.Date.numeric
Loading the Data
load("../data/CAPE_Minder_Rychener_Malsot.RData")
load("../data/NINO34.RData")
load("../data/SRH_Minder_Rychener_Malsot.RData")
Generate PROD
prod = sqrt(cape)*srh
** Create Time Series Objects **
dates <- seq.Date(as.Date("1979-1-1"), as.Date("2015-12-31"), by=1)
feb29ix <- format(as.Date(dates), "%m-%d") == "02-29"
dates <- dates[!feb29ix]
prod_ts = xts(prod, order.by = rep(dates, each=8))
Beginning the Analysis
month_names = c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec")
get_monthly = function(x) {
output = list()
len = nrow(x)
# Get month of element
month = time(x)
month = gsub("....-", "", month)
month = gsub("-..", "", month)
monthlist = unique(month)
for (i in 1:12) {
output[[i]] = x[month == monthlist[i],]
}
names(output) = month_names
return(output)
}
monthly_max = get_monthly(apply.monthly(prod_ts, max))
r = 2
r_monthly_max = get_monthly(apply.monthly(prod_ts, function(x) order(x, decreasing=T)[1:r]))
# get the monthly maxima of prod
m1 = as.data.frame(apply.monthly(prod_ts, max))$V1;
# produce the gumbel qq plot
gumbel_qq = function(x, title="") qqplot(qgumbel(c(1:length(x))/(length(x)+1)),
x,
main = paste("Gumbell Q-Q Plot", title),
xlab = "Theoretical Quantiles",
ylab = "Sample Quantiles") ;
gumbel_qq(m1)
#qqplot(qgumbel(c(1:length(m1))/(length(m1)+1)),m1,main = "Gumbell Q-Q Plot",xlab = "Theoretical Quantiles", ylab = "Sample Quantiles") ;
The qq plot doesn’t fit very well, especially in the lower tail. This is likely due to seasonal dependence.
Fitting GEV to the entire data
# fit gevd with MLE and produce diagnostic plots
fitmax.MLE<-fgev(m1,method="Nelder-Mead")
par(mfrow=c(2,2))
fitmax.MLE
Call: fgev(x = m1, method = "Nelder-Mead")
Deviance: 9272.599
Estimates
loc scale shape
7.519e+03 7.013e+03 4.757e-03
Standard Errors
loc scale shape
421.90201 331.43749 0.05347
Optimization Information
Convergence: successful
Function Evaluations: 171
plot(fitmax.MLE)
Poor fit, probably because the distribution isn’t stationary. This is especially visible in the Probability plot, in which the confidence band is surpassed, indicating a poor fit.
# fit gevd with Bayesian Techniques
# use the previous outputs (rounded) as initial values (use a different shape)
init<-c(1.6e4,4e3,0.1)
# arbitrary priors
mat <- diag(c(10000,10000,100))
pn <- prior.norm(mean=c(0,0,0),cov=mat)
# proposal standard deviation (found by trial and error to get 20%<acceptance rate<40%)
psd<-c(500,0.03,0.02)
# draw 3k samples from markov chain
nit = 30000
thinning = 300
post<-posterior(nit, init=init,prior=pn,lh="gev",data=m1,psd=psd)
# diagnostic plots
MCMC<-mcmc(post[1:nit %% thinning == 0, ], thin=300)
plot(MCMC)
attr(mcmc(post),"ar")
mu sigma xi total
acc.rates 0.24 0.71 0.70 0.55
ext.rates 0.00 0.00 0.01 0.00
#MCMC_stab <- mcmc(post, thin=50, start=1000)
acf(mcmc(post[1:nit %% thinning == 0, ]))
We observe that there seem to be no substantial issues from the autocorrelation plots.
apply(mcmc(post[1:nit %% thinning == 0, ]),2,mean)
mu sigma xi
236.0524535 11181.3806761 -0.1897774
apply(mcmc(post[1:nit %% thinning == 0, ]),2,sd)
mu sigma xi
98.60894002 624.46021143 0.03409059
** Fit with MLE for months separately**
#monthly_fits = lapply(monthly_max,
# function(x) fgev(data.frame(x)[,1], method="Nelder-Mead"))
monthly_fits = list()
error_cases = c(9, 12)
for (i in 1:length(monthly_max)) {
print(i)
if (i %in% error_cases) {
monthly_fits[[i]] = fgev(as.data.frame(monthly_max[[i]])$V1,
method = "Nelder-Mead",
std.err = FALSE)
}
else {
monthly_fits[[i]] = fgev(as.data.frame(monthly_max[[i]])$V1,
method = "Nelder-Mead")
}
}
[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
optimization may not have succeeded
[1] 10
[1] 11
[1] 12
names(monthly_fits) = names(monthly_max)
gumbel_qq(data.frame(monthly_max[[9]])[,1], "September")
gumbel_qq(data.frame(monthly_max[[12]])[,1], "December")
Fit R largest order statistics
largest_10 = get_monthly(apply.monthly(prod_ts, function(x) order(x, decreasing=T)[1:4]))
largest_2_fit = lapply(largest_10, function(x) rlarg.fit(x[ , 1:2]))
$conv
[1] 0
$nllh
[1] 381.7301
$mle
[1] 172.4434267 63.1813324 -0.8246329
$se
[1] 10.3518367 6.2427031 0.1133745
$conv
[1] 0
$nllh
[1] 366.2861
$mle
[1] 155.5451810 51.4464515 -0.7267874
$se
[1] 8.3567636 4.6071290 0.1019313
$conv
[1] 0
$nllh
[1] 348.4217
$mle
[1] 195.3919929 40.8441452 -0.8324715
$se
[1] 6.33288212 4.00545299 0.08820541
$conv
[1] 0
$nllh
[1] 375.033
$mle
[1] 162.9292046 57.8831308 -0.7755148
$se
[1] 9.282767 5.432561 0.099956
$conv
[1] 0
$nllh
[1] 380.4256
$mle
[1] 157.2651459 61.7442705 -0.6554637
$se
[1] 11.5888901 5.2748342 0.1508781
$conv
[1] 0
$nllh
[1] 364.4628
$mle
[1] 104.8006071 51.1390400 -0.1253138
$se
[1] 9.153534 4.541184 0.147210
$conv
[1] 0
$nllh
[1] 386.7889
$mle
[1] 128.6412577 66.3424033 -0.4883028
$se
[1] 11.877915 5.260236 0.139483
$conv
[1] 0
$nllh
[1] 380.7543
$mle
[1] 130.1181620 62.4551520 -0.3960438
$se
[1] 10.9040496 4.7896609 0.1237002
$conv
[1] 0
$nllh
[1] 361.8412
$mle
[1] 111.0849233 49.8496920 -0.2990619
$se
[1] 8.2812262 3.8324125 0.1013576
$conv
[1] 0
$nllh
[1] 377.358
$mle
[1] 75.1622227 56.9409902 0.2522703
$se
[1] 10.0789585 7.2986557 0.1984848
$conv
[1] 0
$nllh
[1] 376.1717
$mle
[1] 106.7510273 60.2707332 -0.1192325
$se
[1] 10.4205289 5.3759827 0.1333882
$conv
[1] 0
$nllh
[1] 377.2483
$mle
[1] 126.4665127 59.0621082 -0.2939312
$se
[1] 9.9836743 4.7019303 0.1170011
lapply(largest_2_fit, rlarg.diag)
Fehler in get(".xts_chob", .plotxtsEnv) :
Objekt '.xts_chob' nicht gefunden
get_se = function(x, ix) {
if (is.null(x$std.err)) 0
else x$std.err[ix]
}
mle_loc = unlist(lapply(monthly_fits, function(x) x$estimate[1]))
mle_loc_se = unlist(lapply(monthly_fits, get_se, 1))
mle_scale = unlist(lapply(monthly_fits, function(x) x$estimate[2]))
mle_scale_se = unlist(lapply(monthly_fits, get_se, 2))
mle_shape = unlist(lapply(monthly_fits, function(x) x$estimate[3]))
mle_shape_se = unlist(lapply(monthly_fits, get_se, 3))
plot_w_err = function(x, y, se, title = NULL) {
max_ix = which.max(y)
min_ix = which.min(y)
plot(x, y,
ylim = c(y[min_ix] - se[min_ix], y[max_ix] + se[max_ix]),
main = title)
arrows(x,y-se,x,y+se, code=3, length=0.02, angle = 90)
}
plot_w_err(1:12, mle_loc, mle_loc_se, "Location Parameter as Estimated with Likelihood")
plot_w_err(1:12, mle_scale, mle_scale_se, "Scale Parameter as Estimated with Likelihood")
plot_w_err(1:12, mle_shape, mle_shape_se, "Shape Parameter as Estimated with Likelihood")
** Fit with Bayesian Methods for months separately**
print(acceptance_rates)
$jan
mu sigma xi total
0.31 0.32 0.33 0.32
$feb
mu sigma xi total
0.24 0.30 0.33 0.29
$mar
mu sigma xi total
0.24 0.32 0.20 0.25
$apr
mu sigma xi total
0.23 0.35 0.28 0.29
$may
mu sigma xi total
0.24 0.29 0.21 0.25
$jun
mu sigma xi total
0.24 0.32 0.22 0.26
$jul
mu sigma xi total
0.23 0.20 0.29 0.24
$aug
mu sigma xi total
0.24 0.21 0.25 0.23
$sep
mu sigma xi total
0.23 0.27 0.20 0.23
$oct
mu sigma xi total
0.24 0.38 0.23 0.28
$nov
mu sigma xi total
0.24 0.34 0.34 0.31
$dec
mu sigma xi total
0.34 0.33 0.30 0.33
TODO-> R largest fit
PART 2 First, we check if the location parameter depends on time using a likelihood ratio test
#monthly_fits = lapply(monthly_max,
# function(x) fgev(data.frame(x)[,1], method="Nelder-Mead"))
ratios = list()
trend = 1:length(as.data.frame(monthly_max[[12]])$V1)
trend = (trend-mean(trend))/sd(trend) # scale and center covariates as recommended
error_cases = c(9, 12)
for (i in 1:length(monthly_max)) {
print(i)
fit_const = fgev(as.data.frame(monthly_max[[i]])$V1,
method = "Nelder-Mead",
std.err = FALSE)
fit_dependant = fgev(as.data.frame(monthly_max[[i]])$V1,
method = "Nelder-Mead",
nsloc = trend,
std.err = FALSE)
ratios[[i]] = fit_const$dev-fit_dependant$dev
}
[1] 1
optimization may not have succeeded
[1] 2
optimization may not have succeeded
[1] 3
[1] 4
optimization may not have succeeded
[1] 5
optimization may not have succeeded
[1] 6
[1] 7
[1] 8
[1] 9
optimization may not have succeeded
[1] 10
[1] 11
optimization may not have succeeded
[1] 12
names(ratios) = names(monthly_max)
chi_95level = qchisq(1-0.05/12,1)
plot(unlist(ratios),main="95% confidence test for time independance, Bonferroni multiple Testing", xlab="Month",ylab="Likelyhood Ratio Statistic")
abline(a=chi_95level,b=0,col="red")
Now, let’s check for independance from ENSO
#monthly_fits = lapply(monthly_max,
# function(x) fgev(data.frame(x)[,1], method="Nelder-Mead"))
ratios = list()
# split nino data into months
n = nino34
dim(n)=c(12,length(as.data.frame(monthly_max[[12]])$V1))
error_cases = c(9, 12)
for (i in 1:length(monthly_max)) {
print(i)
trend = n[i,]
trend = (trend-mean(trend))/sd(trend) # scale and center covariates as recommended
fit_const = fgev(as.data.frame(monthly_max[[i]])$V1,
method = "Nelder-Mead",
std.err = FALSE)
fit_dependant = fgev(as.data.frame(monthly_max[[i]])$V1,
method = "Nelder-Mead",
nsloc = trend,
std.err = FALSE)
ratios[[i]] = fit_const$dev-fit_dependant$dev
}
[1] 1
optimization may not have succeeded
[1] 2
optimization may not have succeeded
[1] 3
[1] 4
optimization may not have succeeded
[1] 5
[1] 6
optimization may not have succeeded
[1] 7
[1] 8
optimization may not have succeeded
[1] 9
optimization may not have succeeded
[1] 10
[1] 11
[1] 12
names(ratios) = names(monthly_max)
chi_95level = qchisq(1-0.05/12,1)
plot(unlist(ratios),main="95% confidence test for independance from ENSO, Bonferroni multiple testing", xlab="Month",ylab="Likelyhood Ratio Statistic")
abline(a=chi_95level,b=0,col="red")
Another method is the chi plot:
# plot the chi plot for dependance analysis
trend = 1:length(as.data.frame(monthly_max[[12]])$V1)
n = nino34
dim(n)=c(12,length(as.data.frame(monthly_max[[12]])$V1))
for (i in 1:length(monthly_max)) {
nino = n[i,]
m_data=as.data.frame(monthly_max[[i]])$V1
dat.m1_month = cbind(m_data,trend);
dat.m1_nino = cbind(m_data,nino);
par(mfrow=c(2,2))
chiplot(dat.m1_month,main1 = "Chi Plot Time",main2 = "Chi Bar Plot Time");
chiplot(dat.m1_nino,main1 = "Chi Plot ENSO",main2 = "Chi Bar Plot ENSO");
}
PART 3 We will now analyse temporal clustering of extremes. For this, we will use the exiplot function from the evd library.
# define a function for getting the extremal indices for each month for a given threshold
monthly_eindex <- function(data, threshold_p, r=0){
ei = list()
for (i in 1:length(data)) {
threshold = quantile(as.data.frame(data[[i]])$V1, threshold_p)
ei[[i]]=exi(as.data.frame(data[[i]])$V1, threshold, r)
}
names(ei) = names(data)
return(ei)
}
ei = monthly_eindex(get_monthly(prod_ts), 0.95)
plot(unlist(ei), main="Extremal Index by Month, 95%-Quantile as Threshold", xlab="Month", ylab="Extremal index")
We observe that the extremal index is 0.25-0.45, we can therefore conclude that we have strong dependance of extremes, with the limiting mean cluster size being roughly from 2 to 4. The clustering has no effect for estimaters based on the (monthly) maximum, but the r largest estimator is influenced by it.
PART 4 First, let’s estimate the 10 year return level using point process approach
monthly_fits_pp = list()
monthly_data = get_monthly(prod_ts)
error_cases = c(5,6,7,8,9,10,11,12)
month_days = c(31,28,31,30,31,30,31,31,30,31,30,31)
for (i in 1:length(monthly_max)) {
print(i)
threshold = quantile(as.data.frame(monthly_data[[i]])$V1, 0.95)
if (i %in% error_cases) {
monthly_fits_pp[[i]] = fpot(as.data.frame(monthly_data[[i]])$V1,
threshold = threshold,
model="pp",
npp = month_days[i]*8,
cmax = TRUE,
r = 1,
std.err = FALSE,
method = "Nelder-Mead")
}
else {
monthly_fits_pp[[i]] = fpot(as.data.frame(monthly_data[[i]])$V1,
threshold = threshold,
model="pp",
npp = month_days[i]*8,
cmax = TRUE,
r = 1,
method = "Nelder-Mead")
}
}
[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
[1] 11
[1] 12
names(monthly_fits_pp) = names(monthly_data)
for(i in 1:12){
par(mfrow=c(2,2))
plot(monthly_fits_pp[[i]])
}
The fit in february has completely failed, and the others are not very good either
We will still estimate the return levels:
return_level = function(x,period=20){
p = 1/period
loc = x$estimate[[1]]
scale = x$estimate[[2]]
shape = x$estimate[[3]]
level = loc + scale*(((-log(1-p))^-shape-1)/shape)
return(level)
}
return_level_20 = lapply(monthly_fits, return_level) # 20 for testing
return_level_100 = lapply(monthly_fits, return_level, period=100)
return_level_50 = lapply(monthly_fits, return_level, period=50)
plot(unlist(return_level_100),main="100 Year Return level, estimated with point process", xlab="Month",ylab="Return Level")
plot(unlist(return_level_50),main="100 Year Return level, estimated with point process", xlab="Month",ylab="Return Level")
TODO: estimat with mcmc Assuming that we have the posterior densities of the markov chains, call theis function to plot return level histograms
return_level_mcmc = function(posterior,period=20){
u = mc.quant(posterior,p=1-1/period,lh="gev")
label_mcmc_rl = sprintf("%s Year Return Level",period)
hist(u,nclass=20,prob=T,xlab=label_mcmc_rl, main = "Return Level Histogram")
}
lapply(monthly_bayes_fit, function(x) return_level_mcmc(x$posterior))
$jan
$breaks
[1] 0 5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 55000 60000 65000
[15] 70000 75000 80000 85000 90000
$counts
[1] 14646 12954 1775 383 131 58 15 17 7 4 3 3 0 1
[15] 1 2 0 1
$density
[1] 9.763675e-05 8.635712e-05 1.183294e-05 2.553248e-06 8.733042e-07 3.866538e-07
[7] 9.999667e-08 1.133296e-07 4.666511e-08 2.666578e-08 1.999933e-08 1.999933e-08
[13] 0.000000e+00 6.666444e-09 6.666444e-09 1.333289e-08 0.000000e+00 6.666444e-09
$mids
[1] 2500 7500 12500 17500 22500 27500 32500 37500 42500 47500 52500 57500 62500 67500
[15] 72500 77500 82500 87500
$xname
[1] "u"
$equidist
[1] TRUE
attr(,"class")
[1] "histogram"
$feb
$breaks
[1] 0 5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 55000 60000 65000
$counts
[1] 1 18276 10221 1129 242 74 32 14 5 2 2 1 2
$density
[1] 6.666444e-09 1.218359e-04 6.813773e-05 7.526416e-06 1.613280e-06 4.933169e-07
[7] 2.133262e-07 9.333022e-08 3.333222e-08 1.333289e-08 1.333289e-08 6.666444e-09
[13] 1.333289e-08
$mids
[1] 2500 7500 12500 17500 22500 27500 32500 37500 42500 47500 52500 57500 62500
$xname
[1] "u"
$equidist
[1] TRUE
attr(,"class")
[1] "histogram"
$mar
$breaks
[1] 4000 6000 8000 10000 12000 14000 16000 18000 20000 22000 24000 26000 28000 30000
[15] 32000 34000 36000 38000 40000 42000 44000 46000 48000 50000 52000 54000
$counts
[1] 1 1 0 0 0 0 103 1680 7830 11503 5998 1980 593 178
[15] 66 35 16 11 2 2 1 0 0 0 1
$density
[1] 1.666611e-08 1.666611e-08 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[7] 1.716609e-06 2.799907e-05 1.304957e-04 1.917103e-04 9.996333e-05 3.299890e-05
[13] 9.883004e-06 2.966568e-06 1.099963e-06 5.833139e-07 2.666578e-07 1.833272e-07
[19] 3.333222e-08 3.333222e-08 1.666611e-08 0.000000e+00 0.000000e+00 0.000000e+00
[25] 1.666611e-08
$mids
[1] 5000 7000 9000 11000 13000 15000 17000 19000 21000 23000 25000 27000 29000 31000
[15] 33000 35000 37000 39000 41000 43000 45000 47000 49000 51000 53000
$xname
[1] "u"
$equidist
[1] TRUE
attr(,"class")
[1] "histogram"
$apr
$breaks
[1] 0 5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 55000 60000 65000
[15] 70000
$counts
[1] 1 3 3 3 88 20707 8892 270 21 7 1 3 1 1
$density
[1] 6.666444e-09 1.999933e-08 1.999933e-08 1.999933e-08 5.866471e-07 1.380421e-04
[7] 5.927802e-05 1.799940e-06 1.399953e-07 4.666511e-08 6.666444e-09 1.999933e-08
[13] 6.666444e-09 6.666444e-09
$mids
[1] 2500 7500 12500 17500 22500 27500 32500 37500 42500 47500 52500 57500 62500 67500
$xname
[1] "u"
$equidist
[1] TRUE
attr(,"class")
[1] "histogram"
$may
$breaks
[1] 0 5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 55000 60000 65000
[15] 70000 75000 80000 85000
$counts
[1] 1 3 3 3 1 1 1979 24216 3640 125 15 7 2 2
[15] 2 0 1
$density
[1] 6.666444e-09 1.999933e-08 1.999933e-08 1.999933e-08 6.666444e-09 6.666444e-09
[7] 1.319289e-05 1.614346e-04 2.426586e-05 8.333056e-07 9.999667e-08 4.666511e-08
[13] 1.333289e-08 1.333289e-08 1.333289e-08 0.000000e+00 6.666444e-09
$mids
[1] 2500 7500 12500 17500 22500 27500 32500 37500 42500 47500 52500 57500 62500 67500
[15] 72500 77500 82500
$xname
[1] "u"
$equidist
[1] TRUE
attr(,"class")
[1] "histogram"
$jun
$breaks
[1] 0 5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 55000 60000 65000
[15] 70000 75000 80000 85000
$counts
[1] 1 3 3 3 1 1 2192 22381 5201 185 17 6 2 2
[15] 2 0 1
$density
[1] 6.666444e-09 1.999933e-08 1.999933e-08 1.999933e-08 6.666444e-09 6.666444e-09
[7] 1.461285e-05 1.492017e-04 3.467218e-05 1.233292e-06 1.133296e-07 3.999867e-08
[13] 1.333289e-08 1.333289e-08 1.333289e-08 0.000000e+00 6.666444e-09
$mids
[1] 2500 7500 12500 17500 22500 27500 32500 37500 42500 47500 52500 57500 62500 67500
[15] 72500 77500 82500
$xname
[1] "u"
$equidist
[1] TRUE
attr(,"class")
[1] "histogram"
$jul
$breaks
[1] 4000 6000 8000 10000 12000 14000 16000 18000 20000 22000 24000 26000 28000 30000
[15] 32000 34000 36000 38000 40000 42000
$counts
[1] 2 3 1 3 2 0 1 0 522 13995 14730 663 54 15
[15] 6 2 1 0 1
$density
[1] 3.333222e-08 4.999833e-08 1.666611e-08 4.999833e-08 3.333222e-08 0.000000e+00
[7] 1.666611e-08 0.000000e+00 8.699710e-06 2.332422e-04 2.454918e-04 1.104963e-05
[13] 8.999700e-07 2.499917e-07 9.999667e-08 3.333222e-08 1.666611e-08 0.000000e+00
[19] 1.666611e-08
$mids
[1] 5000 7000 9000 11000 13000 15000 17000 19000 21000 23000 25000 27000 29000 31000
[15] 33000 35000 37000 39000 41000
$xname
[1] "u"
$equidist
[1] TRUE
attr(,"class")
[1] "histogram"
$aug
$breaks
[1] 4000 6000 8000 10000 12000 14000 16000 18000 20000 22000 24000 26000 28000 30000
[15] 32000 34000 36000 38000 40000 42000 44000
$counts
[1] 1 1 0 0 2 0 0 6 216 3344 12190 10653 2914 518
[15] 115 20 14 2 2 3
$density
[1] 1.666611e-08 1.666611e-08 0.000000e+00 0.000000e+00 3.333222e-08 0.000000e+00
[7] 0.000000e+00 9.999667e-08 3.599880e-06 5.573148e-05 2.031599e-04 1.775441e-04
[13] 4.856505e-05 8.633046e-06 1.916603e-06 3.333222e-07 2.333256e-07 3.333222e-08
[19] 3.333222e-08 4.999833e-08
$mids
[1] 5000 7000 9000 11000 13000 15000 17000 19000 21000 23000 25000 27000 29000 31000
[15] 33000 35000 37000 39000 41000 43000
$xname
[1] "u"
$equidist
[1] TRUE
attr(,"class")
[1] "histogram"
$sep
$breaks
[1] 0 5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 55000 60000 65000
$counts
[1] 1 3 3 15 23255 6680 32 4 0 2 2 3 1
$density
[1] 6.666444e-09 1.999933e-08 1.999933e-08 9.999667e-08 1.550282e-04 4.453185e-05
[7] 2.133262e-07 2.666578e-08 0.000000e+00 1.333289e-08 1.333289e-08 1.999933e-08
[13] 6.666444e-09
$mids
[1] 2500 7500 12500 17500 22500 27500 32500 37500 42500 47500 52500 57500 62500
$xname
[1] "u"
$equidist
[1] TRUE
attr(,"class")
[1] "histogram"
$oct
$breaks
[1] 4000 6000 8000 10000 12000 14000 16000 18000 20000 22000 24000 26000 28000 30000
[15] 32000 34000 36000 38000 40000 42000 44000 46000 48000
$counts
[1] 1 1 0 0 0 30 938 5430 9239 8150 3932 1437 485 204 70 44 20
[18] 9 3 3 1 4
$density
[1] 1.666611e-08 1.666611e-08 0.000000e+00 0.000000e+00 0.000000e+00 4.999833e-07
[7] 1.563281e-05 9.049698e-05 1.539782e-04 1.358288e-04 6.553115e-05 2.394920e-05
[13] 8.083064e-06 3.399887e-06 1.166628e-06 7.333089e-07 3.333222e-07 1.499950e-07
[19] 4.999833e-08 4.999833e-08 1.666611e-08 6.666444e-08
$mids
[1] 5000 7000 9000 11000 13000 15000 17000 19000 21000 23000 25000 27000 29000 31000
[15] 33000 35000 37000 39000 41000 43000 45000 47000
$xname
[1] "u"
$equidist
[1] TRUE
attr(,"class")
[1] "histogram"
$nov
$breaks
[1] 4000 6000 8000 10000 12000 14000 16000 18000 20000 22000 24000 26000 28000 30000
[15] 32000 34000 36000 38000 40000 42000 44000
$counts
[1] 1 102 5217 14247 7181 2090 689 269 101 50 19 14 8 1
[15] 2 5 2 2 0 1
$density
[1] 1.666611e-08 1.699943e-06 8.694710e-05 2.374421e-04 1.196793e-04 3.483217e-05
[7] 1.148295e-05 4.483184e-06 1.683277e-06 8.333056e-07 3.166561e-07 2.333256e-07
[13] 1.333289e-07 1.666611e-08 3.333222e-08 8.333056e-08 3.333222e-08 3.333222e-08
[19] 0.000000e+00 1.666611e-08
$mids
[1] 5000 7000 9000 11000 13000 15000 17000 19000 21000 23000 25000 27000 29000 31000
[15] 33000 35000 37000 39000 41000 43000
$xname
[1] "u"
$equidist
[1] TRUE
attr(,"class")
[1] "histogram"
$dec
$breaks
[1] 0 5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 55000 60000 65000
[15] 70000
$counts
[1] 23735 5642 470 103 16 18 9 1 2 0 0 0 3 2
$density
[1] 1.582281e-04 3.761208e-05 3.133229e-06 6.866438e-07 1.066631e-07 1.199960e-07
[7] 5.999800e-08 6.666444e-09 1.333289e-08 0.000000e+00 0.000000e+00 0.000000e+00
[13] 1.999933e-08 1.333289e-08
$mids
[1] 2500 7500 12500 17500 22500 27500 32500 37500 42500 47500 52500 57500 62500 67500
$xname
[1] "u"
$equidist
[1] TRUE
attr(,"class")
[1] "histogram"